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Abstract 

Let  X :=  L1([0,  ro]),  where  to  represents  the  optical  depth  of  a stel- 
lar atmosphere.  The  weakly  singular  integral  operator  T : X— >X  defined  by 

PV)(r)  = f fj°  Ei(\t  - t'|V(t')  dr', 

where  tu  €]0, 1[  is  the  albedo  of  the  atmosphere  and  E\  denotes  the  first 
exponential-integral  function,  is  such  that  ||T||,  = ro(l  — JE2(to/2)),  where 

E2  denotes  the  second  exponential-integral  function.  If  w is  close  to  1,  and 
To  is  large,  then  ||T||,  is  close  to  1.  In  that  case,  the  transfer  problem 
given  f€X,  find  (p( zX  such  that  Tip  = ip  + / 

is  ill-conditioned,  and  the  convergence  of  the  fixed-point  iteration  ipk+ 1 = Tipk  — f,  which 
is  commonly  used  by  numerical  astronomers,  becomes  prohibitively  slow.  The  purposes 
of  this  work  are  to  approximate  p through  different  sequences  whose  terms  solve  well- 
conditioned  approximate  equations,  and  to  compare  their  efficiency  and  computational 
costs. 


1 Introduction 

For  a given  t0  > 0,  let  g be  a function  defined  on  ]0,  t0]  such  that 


lim  g(r)  = +00, 

t— *0+ 

(1.1) 

5GCo(]0)To])nL1([0,ro]), 

(1.2) 

g(r)  > 0 for  all  r € ]0,to], 

(1.3) 

g is  a decreasing  function  on  ]0,  To]. 

(1.4) 

We  consider  the  integral  operator  T defined  by 

(Tx)(t)  :=  f g{\T-T’\)x{T')dT'. 
Jo 

(1.5) 

rr 0/2 

Theorem  1 T is  a linear  compact  operator  in  L1([0,ro])  and  HT^  = 2 / g(r)  dr. 

Jo 

Proof:  See  [2].  □ 
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For  2 in  the  resolvent  set  of  T,  we  consider  the  Fredholm  equation  of  the  second  kind 

T(p  = z<p  + f.  (1.6) 

Applications  will  concern  the  function  g : ]0,  To]  — ► R given  by 

(17) 

’exp  (-rp) 


w . 


9(t)  ■■=  jEi(t) 


f°°exp< — T Ll) 

where  w e ]0, 1[  and  Ei  is  the  exponential-integral  function  : Ei(t)  :=  / dp, 

J l M 


t > 0.  Ei  is  the  first  function  of  the  sequence  (i5„)„>i,  Ev(t 


>-jf 


5 exp (-tii) 

pv 


dp, 


t > 0,  v > 2,  and  it  is  the  only  one  presenting  a logarithmic  singularity  at  r = 0. 
Following  Theorem  1,  when  g is  defined  by  (1.7),  we  have  HTjl,  = zu[l  — E2(to/2)]  < 1. 

We  recall  that  a bounded  linear  finite  rank  operator  Tn  in  a normed  linear  space  X 
can  be  written  as 


(1.8) 


Tn  > £n,j)Cn,j 

3=1 

where  n € N*,  and,  for  j e [l,nj,  £nj  € A*,  the  topological  adjoint  space  of  X,  and 


°n,3 


ex. 


The  resolution  of  the  approximate  equation 

Tn<Pn  = Zipn  + /,  (1.9) 

where  2 belongs  to  the  resolvent  set  of  Tn,  leads  to  an  n-dimensional  linear  system 

(An  z\n)xn  = bn  (1.10) 

where  l„  is  the  identity  matrix  of  order  n, 

A niilj ) :=  (en,j  J (-■ n,i)i  bn(i)  •=  (/  , I n,i)i  *n{j)  (Vn  i n,j )■  (1-H) 

Once  this  system  is  solved,  the  solution  of  (1.9)  is  given  by 


<Pn  = J !^x„(j)en)j  “ / j ■ 


(1.12) 


We  are  interested  in  refining  approximations  obtained  with  Tn  :=  7 r„T,  where  7r„  is 
a sequence  of  projections  with  finite  rank  n.  A bounded  projection  irn  of  finite  rank  n is 

n 

defined  by  irnx  :=  (x , e*nJ)enj  for  all  x e X,  where  {enj)J=1  is  an  ordered  basis  of 
3= 1 

the  range  of  7r„,  and  (e*  j)"_i  is  an  adjoint  basis  of  the  former  in  X*.  Hence 

n 

Tnx  :=  ^ ',(Tx , enj)en  j , x e X.  (1-13) 

3=1 

We  suppose  that  is  pointwise  convergent  to  the  identity  operator  in  the  Banach  X 
where  the  operator  T is  defined.  Since  T is  compact,  Tn  converges  to  T in  the  operator 
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norm.  Let  R(z)  :=  (I T - zl)  1 be  the  resolvent  of  T at  2.  Then  Rn{z)  (T„  - zl)  1 
exists  for  n large  enough  and  is  uniformly  bounded,  that  is,  there  exists  no  such  that 

c0(z)  :=  sup  ||i?n(z)||  < +00.  (1.14) 

n>«o 

We  develop  an  application  in  the  space  X :=  ^([O,  r0]).  Let  (rnj)"=0  be  a grid  on  [0,  To] 
such  that 


0 — ' Tnj0  <.  Tn>i  • • * < 7n,n— 1 ^ 1 ~n,n  • Ah 


and  set 


hn,j  • — Ti ,7  T~n,j—i  for  j € [1,  • • • ,n]. 
We  define,  for  r G [0,  To], 

„ ._  / 1 ifr€  (Tn!j-i,Tnj) 


p .(T).=  I 1 if -re  (rnii_i,Tn,, 
n,j\  ) ■ | o otherwise 


and,  for  x G Z^QOjTq]), 


(x  > <,j)  '■=  [ x\t ')  dr'. 

hnjJrnj- 1 


The  product  defined  in  (1.18)  is  a special  case  of  the  scalar  product  used  in  equation 
(1.8)  when  a grid  such  as  (1.15)  is  set.  In  this  case  the  operator  in  (1.13)  is  the  operator 
in  (1.8)  if  we  choose  i„:j  = T*e*n  y Let 

Hn  :=  min{/inj-  : j G [1, . . . ,n]},  hn  :=  max{hnJ  : j G [1, . . . , n]},  qn  :=  ^-.(1.19) 

For  quasi-uniform  grids,  there  exists  a constant  q independent  of  n such  that,  for  all  n, 
q<qn-  For  uniform  grids,  qn  = 1 for  all  n. 

Theorem  2 Let  ip  ^ 0 be  the  solution  of  (1-6)  with  T defined  by  (1.5).  Let  tpn  be  the 
solution  of  (1.9)  withTn  defined  by  (1.8)  and  (1.15)-(1.1 7).  Then,  for  n large  enough, 

lly-y-lli  < ggoW  ft  w (1.201 


<5a(£2  fg\T)dr, 

II  i Qn  Jo 

where  co(z)  is  given  by  (1.14)  and  computed  with  the  1 -norm. 

Proof:  See  [2]. 

In  the  case  (1-7),  the  matrix  An  of  the  linear  system  (1.10)  has  entries 
A n{i,j)  ~ TfT—  [ f Ei{\t-  T'\)enj{T')dT'dT, 

and  the  second  member  bn  has  entries 

777  A1"'*  fT° 

bn(*):=xr— / / E1(\T-T,\)f(T,)dT,dT. 

zrin,i  yTni_jJo 
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For  more  details,  see  [3].  An  application  to  the  transfer  problem  in  astrophysics  gives 
(1.6)  with  z = 1,  and  as  free  term, 


f(T\  ._  / -1  if  0 < r < ro/2, 
0 if  To 


' To/ 2 < T < To, 


(1.23) 


which  describes  a sudden  drop  of  the  temperature  on  the  r = 7o/2  layer  of  the  atmo- 
sphere. For  further  details  on  the  physical  model,  see  [4]. 

2 Iterative  refinement  of  approximate  solutions 

To  attain  a given  precision  on  the  approximate  solution  </?„,  it  may  be  necessary  that  the 
largest  grid  step  hn  be  so  small  that  the  dimension  of  the  corresponding  linear  system 
will  be  prohibitively  large  from  a computational  point  of  view.  Not  only  the  algorithm’s 
stability  becomes  poor  but  also  the  condition  number  of  the  matrix  may  increase  if  its 
size  increases.  Refinement  schemes  allow  us  to  attain  iteratively  the  exact  solution  of  a 
large  scale  linear  system  by  means  of  the  resolution  of  a sequence  of  linear  systems  of 
moderate  fixed  size.  Let  us  consider  the  general  framework  of  a complex  Banach  space 
X and  a linear  compact  operator  T : X — > X.  If  2 is  in  the  resolvent  set  of  T,  then 
2^0.  Let  Tn  be  a sequence  of  linear  bounded  operators  in  X such  that  \\T  — Tn\\  — > 0 
in  the  operator  norm.  Then,  for  n large  enough,  2 belongs  to  the  resolvent  set  of  Tn  and 
Rn(z)  is  norm-convergent  to  R(z). 

The  most  elementary  way  to  refine  the  approximate  solution  <pn  :=  Rn(z)f  is  the 
following. 


Scheme  A 


X{0)  •■=  <Pn, 


x(k+l)  x(k)  _ Rn(z)(Tx(k) 


zx 


(k)  _ 


/),  k>  0. 


(2.1) 


We  can  interpret  Rn(z)  as  an  approximation  of  the  inverse  of  the  Frechet  derivative  of 
the  affine  operator  x (T  — zl)x  - f,  the  exact  one  being  R(z).  Since  R(z)  satisfies  the 
identities 

R(z)  = -z(R(z)T-I)=-z(TR(z)-I)  (2.2) 

two  new  different  approximations  of  R(z)  are  thus  motivated, 

Rn(z)  :=  ~(Rn(z)T  - I),  X{z)  :=  -(! TRn(z ) - I).  (2.3) 

2 2 

These  approximate  resolvent  operators  lead  to  the  following  iterative  refinement  schemes, 


Scheme  B 


Scheme  C 


Z(°> 

■=  Rn(z)f, 

x(k+l) 

:=  - Rn(z){Tx^  - 

- ZX^  - /), 

k > 0, 

&°> 

■=  Rn(z)f, 

x(k+l) 

:=  x(k)  - Rn(z)(TxW  - 

-zx^-f), 

k > 0. 

(2.4) 


(2.5) 
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Since  the  computation  of  residuals  which  tend  to  zero,  as  well  as  the  resolution  of  almost 
homogeneous  linear  systems  may  be  unstable,  the  following  theorems  are  interesting  for 
algorithmic  purposes. 

Theorem  3 In  (2.1),  = s<°>  + Rn{z)(Tn  - T)x for  k > 0. 

Theorem  4 In  (2.4),  £(fc+1)  = £(0)  + -Rn(z){Tn  - T)Tx<*>  for  fc  > 0. 

z 

Theorem  5 In  (2.5),  x ^ = P + -TRn(z)(Tn  - T)x for  fc  > 0. 

z 

Proof:  For  each  fc  > 0,  in  (3), 

X(k+1)  _ x(k)  _ Rn(z)(TxW  - zxW  - /) 

- xM -Rn(z){T-Tn  + Tn-zI)xW + xW 
= xM  + Rn(z)(Tn-T)xW 


For  (4)  and  (5),  the  proof  follows  the  same  idea  but  it  is  technically  more  complicated. 

□ 

In  our  application  to  the  transfer  equation  in  astrophysics,  T is  defined  by  (1.5)  with 
g given  by  (1.7),  and  the  equation  (1.6)  has  z = 1. 

3 Numerical  computations 

The  iterative  refinement  schemes  allow  us  to  obtain  the  exact  solution  of  a large  scale 
linear  system  by  solving  a sequence  of  moderate  fixed  size  ones.  Each  of  the  three  iterative 
refinement  schemes  presented  in  this  work  are  based  on  an  approximation,  say  Gn(z), 
of  the  resolvent  operator  R(z).  Their  common  structure  is  the  following. 

€(0)  “ Gn(z)f,  ( , 

^k+i)  £(o  ) + (j_Gn(*)(T  _*/)){(*>,  fc  > 0. 

Theorem  6 Letci(z)  8co(z)max{l,  ||T,||1/|jz:|},  and  {^)k>o  be  any  of  the  sequences 
(2.1),  (2.4)  or  (2.5).  Then 

U(k)~Tl 


< (flM  t’ g(T)dr)k*' , i>0. 
Ml  ~ Wn  Jo  W / 


Proof:  Let  us  prove  the  bound  for  the  sequence  defined  by  (2.1).  For  the  other  two, 
the  arguments  are  similar.  Using  Theorem  3,  we  have 

xW-  <p  = {Rn(z)(Tn  - T))V0)  - <P), 

= Rn(z)(T-Tn)g>. 

Hence, 


r(*) 


^||1<||(JRn(z)(T-Tn))fc+1||1M|1, 


and,  in  [2],  we  have  shown  that  \\Rn(z)(Tn  — T)||,  < 


P 9(T)iT. 

Jo 


□ 


Qn 


All  the  schemes  need  evaluations  of  T at  some  prescribed  functions  of  X.  In  practice 
T is  not  used  for  this  purpose  but  an  operator  Tm  of  the  sequence  (T„)i,>i  is  used  instead, 
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where  m > n.  We  consider  the  kernel  g defined  by  (1.7)  and  the  free  term  / defined 
by  (1.23).  Table  1 gives  the  number  of  iterations  performed  by  each  scheme  for  several 
values  of  w in  order  to  obtain  a first  relative  residual  less  than  or  equal  to  10~12,  when 
a quasi-uniform  grid  (r„,j)£_ 0 is  built  such  that  u is  a multiple  of  10,  r0  = 1000, 

f £ « <e[i,...,|].  ■ 


n = 200,  m = 1000,  and  hvi  := 


£ “ + l f], 


19.  _i_  1 Bill 

2j,  11  * t [2  + ■*■>•••  > 10 J ’ 


4 t0  . 


if  i € [%  + 


Albedo 

Scheme  A 
(2.1) 

Scheme  B 
(2.4) 

Scheme  C 
(2.5) 

0.750 

29 

15 

14 

0.990 

46 

27 

26 

0.999 

385 

196 

195 

TAB  1.  Number  of  iterations. 


Figures  1,  2 and  3 show  the  last  iterate  of  all  schemes,  as  well  as  the  corresponding 
convergence  histories,  for  w e {0.750,0.990,0.999}.  As  we  can  see,  the  schemes  B and 
C are  much  faster  than  Atkinson’s  formula  A,  specially  when  the  albedo  is  close  to  1.  In 
the  latter  situation  a wider  boundary  layer  arises  at  the  left  of  the  atmosphere,  and  the 
decay  at  the  middle  point  takes  place  along  a wider  subinterval. 

A survey  on  different  discretization  methods  for  integral  operators  can  be  found  in 
[1],  with  special  emphasis  on  spectral  applications.  In  what  concerns  condition  number 
of  associated  linear  systems,  the  reader  is  refered  to  [7],  [5]  and  [6]. 
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Solution  Residual 


Fig.  1.  Solution  and  convergence  history  for  zo  — 0.750:  Scheme  A — dashed  line, 
Scheme  B — dotted  line,  Scheme  C — solid  line. 


Solution  Residual 


Fig.  2.  Solution  and  convergence  history  for  w = 0.990:  Scheme  A — dashed  line, 
Scheme  B — dotted  line,  Scheme  C — solid  line. 
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Solution  Residual 


Fig.  3.  Solution  and  convergence  history  for  w = 0.999:  Scheme  A — dashed  line, 
Scheme  B — dotted  line,  Scheme  C — solid  line. 


